#include "RectangleReceiverRectangleGridRayTracing.cuh"
#include <algorithm/rayTracing/dataLoader/CloudLoader.h>
#include <algorithm/rayTracing/rayTracingDataStructure/receiver/CudaRectangleReceiver.cuh>
#include <algorithm/rayTracing/rayTracingPipeline/quasiMonteCarloRayTracing/Grid/rectGridDDA.cuh>
#include <algorithm/rayTracing/rayTracingPipeline/quasiMonteCarloRayTracing/Receiver/rectangleReceiverIntersection.cuh>
#include <algorithm/shared/tool/randomNumberGenerator/RandomGenerator.cuh>

// 考虑三次法相扰动

using namespace solar;
using namespace solar::cuda;

__global__ static void
map_raytracing(SunrayArgument sunrayArgument, CudaRectangleReceiver rectangleReceiver,
               CudaRectangleGrid rectangleGrid, HeliostatArgument heliostatArgument,
               float3* d_subHeliostat_vertexe, float factor, bool res, CloudLoader cloud)
{
    int myId = global_func::getThreadID(); // 线程ID
    int numberOfMicroHeliostats =
        heliostatArgument.micronums * heliostatArgument.helionums; // 所有微定日镜个数
    if (myId >= numberOfMicroHeliostats * sunrayArgument.numberOfLightsPerGroup)
    {
        return; // id比光线数目大，返回
    }
    int microsize = heliostatArgument.micronums;                       // 每个定日镜的微定日镜个数
    int no = myId / microsize / sunrayArgument.numberOfLightsPerGroup; // 第几号镜子

    // int is_total_reflect=0;
    float depth_ = 0.003f;
    float refract_ = 1.523f;
    int helio_no =
        myId / sunrayArgument.numberOfLightsPerGroup / heliostatArgument.micronums; // 第几个定日镜
    int micro_no = myId / sunrayArgument.numberOfLightsPerGroup %
                   heliostatArgument.micronums; // 该定日镜的第几个微定日镜

    int substart = helio_no * 3;

    float3 helio_pos =
        (d_subHeliostat_vertexe[substart] + d_subHeliostat_vertexe[substart + 2]) / 2;
    float3 helio_norm = normalize(
        cross(d_subHeliostat_vertexe[substart] - d_subHeliostat_vertexe[substart + 1],
              d_subHeliostat_vertexe[substart + 2] - d_subHeliostat_vertexe[substart + 1]));

    int col_no = micro_no / heliostatArgument.col;
    int row_no = micro_no % heliostatArgument.col;
    float3 micropos = helio_pos;
    float3 micronorm = helio_norm;

    micropos.z = col_no * heliostatArgument.pixel + heliostatArgument.pixel / 2 -
                 heliostatArgument.sizex / 2;
    micropos.y = heliostatArgument.sizey / 2;
    micropos.x = row_no * heliostatArgument.pixel + heliostatArgument.pixel / 2 -
                 heliostatArgument.sizez / 2;

    float3 local = micropos;
    local = global_func::local2world_rotate(local, micronorm);
    local = global_func::translate(local, helio_pos);
    micropos = local;

    int address = myId % sunrayArgument.pool_size;

    // 采样数组的采样值
    float3 dir = global_func::local2world_rotate(sunrayArgument.d_samplelights[address],
                                                 -sunrayArgument.sunray_direction);
    // 入射光方向（世界坐标系）
    // dir = -sunrayArgument.sunray_direction;
    // 计算第helio_no个定日镜的第micro_no个微定日镜的法向和坐标

    float3 origin = micropos; // 微定日镜中心坐标

    if (rectGridDDA::collision_(origin, dir, rectangleGrid, d_subHeliostat_vertexe,
                                heliostatArgument, no))
    {
        return;
    }
    float3 normal = micronorm; // 微定日镜法向方向，局部坐标系
    //    address = (myId % sunrayArgument.numberOfLightsPerGroup) % sunrayArgument.pool_size;
    float3 normal1 = global_func::local2world_rotate(sunrayArgument.d_perturbations[address],
                                                     normal); // 微定日镜法向方向，世界坐标系
    // normal1 = normal;//微定日镜法向方向，世界坐标系

    normal1 = normalize(normal1);
    dir = normalize(dir); // 入射角dir.y>0入社向量的逆方向(&rectangleReceiver)->

    float theta; // 入射光与法向之间的夹角
    if (std::abs(dot(dir, normal1) - 1.0f) < Epsilon)
        theta = 0.0F;
    else
        theta = acos(dot(dir, normal1));

    // 先和其他定日镜求交,看是否被遮挡
    // 如果没有被遮挡,求光线与云的交点,如果在数组中,取数组中的值,不在,取1.返回值是float值
    float3 dir_n = make_float3(0.0, 1.0, 0.0);
    float cos_theta = dot(dir, dir_n);
    float len = cloud.start_point.y / cos_theta;
    float3 dirs = dir * len;
    float2 M = make_float2(origin.x + dirs.x, origin.z + dirs.z);
    float times = cloud.start_point.y / -sunrayArgument.sunray_direction.y;
    float2 cloudcenter = make_float2(-sunrayArgument.sunray_direction.x * times,
                                     -sunrayArgument.sunray_direction.z * times);
    float2 cloudstart =
        make_float2(cloudcenter.x + cloud.start_point.x, cloudcenter.y + cloud.start_point.z);
    float my = M.x - cloudstart.x; // M.x, M.y 是M点在世界坐标系的（x,z)坐标值
    float mx = M.y - cloudstart.y; // 需要转换成 cloud map局部坐标
    /*世界坐标  --> X
     *         |
     *         |
     *         Z
     * cloud map 局部坐标  --> y
     *                    |
     *                    |
     *                    x
     * M.x-cloud.start_point.x M点在局部坐标系的第几列
     * M.y-cloud.start_point.z M点在局部坐标系的第几行
     * A    Boat3 &dir, float &t, float &u, int &panel_no,int &tube_no,int &angle_no,int
     &ifinterspace,int &interspace_no){
        //If the origin is inside the cylinder, it won't intersect with it
     *    M
     * D    C
     * 根据A,B,C,D四点的值双线性插值得到M的值
     */
    int Ax = floor(mx / cloud.pixel); // A在cloud map中第几行
    int Ay = floor(my / cloud.pixel); // A在cloud map的第几列
    float decay = 1.0;
    if (Ax <= 0 || Ax > cloud.x_size || Ay < 0 || Ay > cloud.y_size)
    {
        decay = 1.0;
        //      printf("times %f cloudstart %f %f %f %f \n",times,cloudstart.x,cloudstart.y,mx,my);
    }
    else
    {
        //        printf("my %f mx %f\n",my,mx);
        /*世界坐标  --> X
         *         |
         *         |
         *         Z
         * cloud map 局部坐标  --> y
         *                    |
         *                    |
         *                    x
         * M.x-cloud.start_point.x M点在局部坐标系的第几列
         * M.y-cloud.start_point.z M点在局部坐标系的第几行
         * A    B
         *    M
         * D
         * 根据A,B,C,D四点的值双线性插值得到M的值
         */

        int A_pos = Ay + Ax * cloud.x_size;

        float Q11 = cloud.d_cloudMap[A_pos];                    // A
        float Q12 = cloud.d_cloudMap[A_pos + 1];                // B
        float Q21 = cloud.d_cloudMap[A_pos + cloud.x_size];     // D
        float Q22 = cloud.d_cloudMap[A_pos + cloud.x_size + 1]; // C
        // A x1,y1
        // B x1,y2
        // C x2,y2
        // D x2,y1
        float x1 = Ax * cloud.pixel;
        float x2 = (Ax + 1) * cloud.pixel;
        float y1 = Ay * cloud.pixel;
        float y2 = (Ay + 1) * cloud.pixel;

        decay = (Q11 * (x2 - mx) * (y2 - my) + Q21 * (mx - x1) * (y2 - my) +
                 Q12 * (x2 - mx) * (my - y1) + Q22 * (mx - x1) * (my - y1)) /
                cloud.pixel / cloud.pixel;
        decay = 1.0 - decay;
    }

    // bool res=1;//0是反射，1是折射
    if (res)
    {
        // 夹角小于全内反射的临界角，折射,调研得到玻璃的全内反射临界角大约为PI/5

        float3 dir_1 = -dir;
        float eta = 1 / refract_;
        float k = 1.0 - eta * eta * (1.0f - dot(normal1, dir_1) * dot(normal1, dir_1));

        float theta_2 = asin(sin(theta) / refract_); // 折射角
        float3 dir2 = normalize(eta * dir_1 - (eta * dot(normal1, dir_1) + sqrt(k)) * normal1);

        float3 dir22 = dir2 * depth_ / cos(theta_2);
        int ran = int(sunrayArgument.pool_size / 3);
        // printf("sunrayArgument.numberOfLightsPerGroup=%d\n",sunrayArgument.numberOfLightsPerGroup);
        // 计算第一次折射向量的反射
        float3 normal2 = normalize(global_func::local2world_rotate(
            sunrayArgument.d_perturbations[(address + ran) % sunrayArgument.pool_size], normal));
        // normal2=normalize(normal);

        float3 dir3 = normalize(reflect(dir2, normal2));
        float3 dir33 = dir3 * depth_ / cos(theta_2);

        // 折射向量的反射vector
        float3 origin1 = origin + dir22 + dir33;

        // 实际反射的原点
        float3 normal3 = normalize(global_func::local2world_rotate(
            sunrayArgument.d_perturbations[(address + ran * 2) % sunrayArgument.pool_size],
            normal));
        // normal3=normalize(normal);

        float theta_refract; // 入射光与法向之间的夹角
        if (abs(dot(dir3, normal3) - 1.0f) < Epsilon)
        {
            theta_refract = 0.0f;
        }
        else
        {
            theta_refract = std::acos(dot(dir3, normal3));
        }
        float theta_reflect = std::asin(1.0F / refract_);

        if (theta_refract < theta_reflect)
        {
            // 反射向量与玻璃面反面法线夹角小于全内反射角
            float eta_1 = refract_;
            float k_1 = 1.0 - eta_1 * eta_1 * (1.0f - dot(-normal3, dir3) * dot(-normal3, dir3));

            float3 dir4 =
                normalize(eta_1 * dir3 - (eta_1 * dot(-normal3, dir3) + sqrt(k_1)) * (-normal3));
            // 折射向量
            dir = dir4;
            origin = origin1;
            normal1 = normal3;
        }
        else
        {
            return;
        }
    }
    else
    {
        dir = normalize(reflect(-dir, normal1));
        normal = normal1;
    }
    if (rectGridDDA::collision_(origin, dir, rectangleGrid, d_subHeliostat_vertexe,
                                heliostatArgument, no))
    {
        return;
    }
    float t;
    float u;
    float v;

    if (!rectangleReceiver.GIntersect(origin, dir, t, u, v))
    {
        return;
    }

    // Step 2: Calculate the energy of the light
    float eta = 0.0;
    if (t <= 1000.0F)
    {
        eta = 0.99331F - 0.0001176F * t + 1.97F * (1e-8F) * t * t;
    }
    else
    {
        eta = expf(-0.0001106F * t);
    }
    float energy = fabsf(dot(dir, normal1)) * eta * factor * decay;

    int nx = rectangleReceiver.resolution_.x;
    int ny = rectangleReceiver.resolution_.y;
    float* imagess = rectangleReceiver.getDeviceImage();
    int2 row_col = make_int2(u * ny, v * nx);
    int addresses = row_col.x * nx + row_col.y; // col_row.y + col_row.x * resolution.y
    atomicAdd(&(imagess[addresses]), energy);   // CUDA atomic
}

void solar::cuda::RectangleReceiverRectangleGridRayTracing(
    SunrayArgument* sunrayArgument, CudaRectangleReceiver* rectangleReceiver,
    CudaRectangleGrid* rectGrid, HeliostatArgument* heliostatArgument,
    float3* d_subHeliostat_vertexe, float factor, bool res, CloudLoader* cloud, cudaStream_t stream)
{
    int n_threads = 512;
    dim3 n_blocks;
    /* //float mean_ =0.0022f;
     //float err_ = 0.0015f;

     //float *gaussian = new float[heliostat_nums];

     //Step 3: Generate theta and phi
     //RandomGenerator::cpuGaussian(gaussian, mean_, err_, heliostat_nums);

     //float *d_phi= nullptr;
     //checkCudaErrors(cudaMalloc((void **) &d_phi,  sizeof(float) * heliostat_nums));

     //RandomGenerator::gpuUniform(d_phi,heliostat_nums);


     //float *d_gaussian= nullptr;
     //checkCudaErrors(cudaMalloc((void **) &d_gaussian,  sizeof(float) * heliostat_nums));
     //global_func::cpu2gpu(d_gaussian,gaussian,heliostat_nums);
 */
    global_func::setThreadBlocks(n_blocks, n_threads,
                                 heliostatArgument->helionums * heliostatArgument->micronums *
                                     sunrayArgument->numberOfLightsPerGroup,
                                 true);

    // printf("num of rays= %lld \n", heliostatArgument->helionums * heliostatArgument->micronums *
    //                                    sunrayArgument->numberOfLightsPerGroup);
    // printf("nBlocks = %d %d %d, nThreads = %d\n", n_blocks.x, n_blocks.y, n_blocks.z, n_threads);

    map_raytracing<<<n_blocks, n_threads, 0, stream>>>(*sunrayArgument, *rectangleReceiver, *rectGrid,
                                            *heliostatArgument, d_subHeliostat_vertexe, factor, res,
                                            *cloud);
    // cudaDeviceSynchronize();
    checkCudaErrors(cudaGetLastError());
}
